MBS COVID
Caroline X.Gao, Centre for Youth Mental Health, University of Melbourne, Melbourne, Victoria, Austrlia; Orygen, Parkview, Victoria, Austrlia; School of Public Health and Preventive Medicine, Monash University, Melbourne, Victoria, Austrlia
December 20, 2021
#load libraries
library(pacman)
p_load("tidyverse","ggpubr", "ggspatial", "cowplot", "here","gridExtra","rgeos","rgdal","rmapshaper", "MASS", "nlme", "mgcv","data.table", "kableExtra","gtsummary","Hmisc","arsenal","imputeTS","zoo") source(here::here("COVID","print_table_flextable.R"))Introduction
This file includes all R code for the analysis carried out for the paper “Impact of COVID-19 pandemic on mental health medical service utilization in Australian youth: an interrupted time series analysis”.
Item table
Items<- read_csv(here::here("Data", "Item code.csv")) %>%
mutate(`Group Description`=str_to_sentence(`Group Description`)) %>%
dplyr::select(`MBS Group`,`Group Description`,`MBS Requested Item`) %>%
group_by(`MBS Group`,`Group Description`) %>%
mutate(Items = paste0(`MBS Requested Item`, collapse = ", ")) %>%
dplyr::select( -`MBS Requested Item`) %>%
unique() %>%
mutate(`Group Description`=str_replace(`Group Description`,"Gp", "GP")) %>%
mutate(`Group Description`=str_replace(`Group Description`,"gp", "GP")) %>%
mutate(`Group Description`=str_replace(`Group Description`,"Covid", "COVID"))
print_table(Items, align = c("l", "l", "l"), cap = "MBS items included", capname = "Items", colwidth = c(1,5,10))MBS Group | Group Description | Items |
A06 | Group therapy | 170, 171, 172 |
A07 | Acupuncture and non-specialist practitioner items | 221, 222, 223, 272, 276, 277, 279, 281, 282, 283, 285, 286, 287, 371, 372, 894, 896, 898 |
A08 | Consultant psychiatrist attendances to which no other item applies | 288, 289, 291, 293, 296, 297, 299, 300, 302, 304, 306, 308, 310, 312, 314, 316, 318, 319, 320, 322, 324, 326, 328, 330, 332, 334, 336, 338, 342, 344, 346, 348, 350, 352, 353, 355, 356, 357, 358, 359, 361, 364, 366, 367, 369, 370 |
A15 | GP management plans, team care arrangements, multidisciplinary care plans | 855, 857, 858, 861, 864, 866 |
A30 | Medical practitioner (including a general practitioner, specialist or consultant physician) telehealth attendances | 2121, 2150, 2196 |
A20 | GP mental health treatment | 2700, 2701, 2702, 2710, 2712, 2713, 2715, 2717, 2719, 2721, 2723, 2725, 2727, 2729, 2731 |
M03 | Allied health services | 10956, 10968 |
T01 | Miscellaneous therapeutic procedures | 14224 |
M06 | Psychological therapy services | 80000, 80001, 80005, 80010, 80011, 80015, 80020, 80021 |
M07 | Focussed psychological strategies (allied mental health) | 80100, 80101, 80105, 80110, 80111, 80115, 80120, 80121, 80125, 80126, 80130, 80135, 80136, 80140, 80145, 80146, 80150, 80151, 80155, 80160, 80161, 80165, 80170, 80171 |
M11 | Allied health services for indigenous australians who have had a health check | 81325, 81355 |
M10 | Autism, pervasive developmental disorder and disability services | 82000, 82005, 82010, 82015, 82020, 82025, 82030, 82035 |
M16 | Eating disorders services | 82350, 82351, 82352, 82353, 82354, 82355, 82356, 82357, 82358, 82359, 82360, 82361, 82362, 82363, 82364, 82365, 82366, 82367, 82368, 82369, 82370, 82371, 82372, 82373, 82374, 82375, 82376, 82377, 82378, 82379, 82380, 82381, 82382, 82383 |
A36 | Eating disorder services | 90250, 90251, 90252, 90253, 90254, 90255, 90256, 90257, 90260, 90261, 90262, 90263, 90264, 90265, 90266, 90267, 90268, 90269, 90271, 90272, 90273, 90274, 90275, 90276, 90277, 90278, 90279, 90280, 90281, 90282 |
M17 | Bushfire recovery access initiative - psychologist services and allied health focussed psychological strategies | 91000, 91001, 91005, 91010, 91011, 91015, 91100, 91101, 91105, 91110, 91111, 91115, 91125, 91126, 91130, 91135, 91136, 91140, 91150, 91151, 91155, 91160, 91161, 91165 |
M18 | COVID-19 allied health telehealth services | 91166, 91167, 91169, 91170, 91172, 91173, 91175, 91176, 91181, 91182, 91183, 91184, 91185, 91186, 91187, 91188, 93032, 93033, 93035, 93036, 93040, 93041, 93043, 93044, 93074, 93076, 93079, 93084, 93087, 93092, 93095, 93100, 93103, 93108, 93110, 93113, 93118, 93121, 93126, 93129, 93134, 93137 |
A39 | Bushfire recovery access initiative - GP and medical practitioner focussed psychological strategies | 91283, 91285, 91286, 91287, 91371, 91372, 91721, 91723, 91725, 91727, 91729, 91731 |
A40 | COVID-19 services | 91818, 91819, 91820, 91821, 91827, 91828, 91829, 91830, 91831, 91837, 91838, 91839, 91840, 91841, 91842, 91843, 91844, 91845, 92112, 92113, 92114, 92115, 92116, 92117, 92118, 92119, 92120, 92121, 92122, 92123, 92124, 92125, 92126, 92127, 92128, 92129, 92130, 92131, 92132, 92133, 92134, 92135, 92146, 92147, 92148, 92149, 92150, 92151, 92152, 92153, 92154, 92155, 92156, 92157, 92158, 92159, 92160, 92161, 92162, 92163, 92166, 92167, 92170, 92171, 92172, 92173, 92176, 92177, 92178, 92179, 92182, 92184, 92186, 92188, 92194, 92196, 92198, 92200, 92434, 92435, 92436, 92437, 92455, 92456, 92457, 92458, 92459, 92460, 92474, 92475, 92476, 92477, 92495, 92496, 92497, 92498, 92499, 92500 |
A41 | COVID-19 additional focussed psychological strategies | 93300, 93301, 93302, 93303, 93304, 93305, 93306, 93307, 93308, 93309, 93310, 93311 |
M25 | COVID-19 additional psychological therapy services | 93330, 93331, 93332, 93333, 93334, 93335 |
M26 | COVID-19 additional focussed psychological strategies (allied mental health) | 93350, 93351, 93352, 93353, 93354, 93355, 93356, 93357, 93358, 93359, 93360, 93361, 93362, 93363, 93364, 93365, 93366, 93367 |
Import data
The data include service uptake count data by SA3, year and population subgroups. Other avaliable inforamtion include SA3 region type (regional or metropolitan), state (Victoria vs. other states), estimated residential population (ERP) and the 2016 Socio-economic Indexes for Areas (SEIFA) scores for each SA3 (estimated based on population weighted average of SA2 scores) were used in this analysis. The measure of SEIFA used was the Index of Relative Socio-economic Advantage and Disadvantage (IRSAD).
dta <- readRDS( here("COVID","COVID_cleaned.rds"))
map <- readRDS(here("Data", "Shapefile_SA3.rds"))MBS service rates by SA3
This heat map displays the annual number of services accessed per 1000 people for the combined youth population in 200 for each SA3 in Australia. Separate figures show the heat maps specifically for Greater Brisbane, Greater Sydney and Greater Melbourne.
library(raster)
# state boundary
State <- raster::aggregate(map, "STE_NAME16")
# calculate annual service access rates
rates <- dta %>%
filter( Population == "All young people") %>%
group_by(SA3_CODE16, Year) %>%
summarise(Services = sum(Services), ERP = mean(ERP)) %>%
mutate(Rates = Services * 1000 / ERP) %>%
mutate(SA3_CODE16 = as.character(SA3_CODE16))map@data <- map@data %>%
left_join(rates[rates$Year == 2020, ], by = "SA3_CODE16")
Aus <- ggplot() +
# plot pologon
ggspatial::annotation_spatial(map) +
layer_spatial(map, aes(fill = Rates), colour = "gray68") +
layer_spatial(State, colour = "gray20", fill = NA) +
annotation_scale(location = "bl", width_hint = 0.5) +
annotation_north_arrow(
location = "bl", which_north = "true",
pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering
) +
scale_fill_continuous(low = "white", high = "#e34a33", limits = c(0, 1179)) +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(0.15, 0.85),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.border = element_blank(),
) +
labs(fill = "Number of services \n per 1000 population") +
coord_sf(xlim = c(109, 165), ylim = c(-9, -45))
Mel <- ggplot() +
# plot pologon
ggspatial::annotation_spatial(map) +
layer_spatial(map, aes(fill = Rates), colour = "gray68") +
scale_fill_continuous(low = "white", high = "#e34a33", limits = c(0, 1179)) +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
) +
coord_sf(xlim = c(144.45, 145.55), ylim = c(-37.5, -38.3))
Syd <- ggplot() +
# plot pologon
ggspatial::annotation_spatial(map) +
layer_spatial(map, aes(fill = Rates), colour = "gray68") +
scale_fill_continuous(low = "white", high = "#e34a33", limits = c(0, 1179)) +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
) +
coord_sf(xlim = c(150.67, 151.32), ylim = c(-33.65, -34.15))
Bri <- ggplot() +
# plot pologon
ggspatial::annotation_spatial(map) +
layer_spatial(map, aes(fill = Rates), colour = "gray68") +
scale_fill_continuous(low = "white", high = "#e34a33", limits = c(0, 1179)) +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
) +
coord_sf(xlim = c(152.65, 153.5), ylim = c(-27.05, -27.75))
arrows <- data.frame(xs = c(0.62, 0.71, 0.74), xe = c(0.79, 0.79, 0.79), ys = c(0.23, 0.33, 0.5), ye = c(0.18, 0.43, 0.67))
ggdraw(Aus) +
# add AUD map
draw_plot(Bri, x = 0.75, y = 0.55, width = 0.25, height = 0.25) +
draw_plot(Syd, x = 0.75, y = 0.3, width = 0.25, height = 0.25) +
draw_plot(Mel, x = 0.75, y = 0.05, width = 0.25, height = 0.25) +
geom_segment(aes(x = xs, y = ys, xend = xe, yend = ye),
data = arrows,
arrow = arrow(), lineend = "round"
)
Imputation
Series were removed from the dataset if they had fewer than 15 non-missing data points. The remaining missing values were imputed using linear interpolation, generating integers between 0 and 10.
# extract missing
dta <- dta %>%
group_by(SA3_CODE16, Population) %>%
mutate(missing = sum(is.na(Services))) %>%
ungroup()
# check missing
# table(dta$missing)
set.seed(12345)
dta_imputed<-dta %>%
#keep if the series had less than 15 missing data points
filter( missing <= 15) %>%
# Remove few data points for newly established areas where ERP was too small
filter(ERP>5)%>%
# fill missing with random number between 0 and 9
mutate(Services=na_random(Services,lower_bound= -0.5, upper_bound=9.5)) %>%
mutate(Services=round(Services))
#prepare data for GLM
dta_imputed <- dta_imputed %>%
mutate(YearQuarter=as.numeric(YearQuarter),
COVID = as.numeric(Year==2020),
ID = paste0(SA3_CODE16,"_",Population)) %>%
arrange(ID,YearQuarter)
# checking
# hist(dta$Services[dta$Services<50])
# hist(dta_imputed$Services[dta_imputed$Services<50])
# length(unique(dta_imputed$SA3_CODE16))Line plot
The line plot below presents the quarterly number of services delivered per 1000 people for each SA3 by population subgroup, along with the mean rate at each quarter for that subpopulation.
dta %>%
mutate(YearQuarter = as.Date(YearQuarter)) %>%
group_by(YearQuarter,Population) %>%
mutate(Mean = mean(Rates,na.rm = TRUE) ) %>%
ungroup() %>%
ggplot(aes(x = YearQuarter, y = Rates, group = SA3_CODE16)) +
geom_line(alpha = 0.5, col = "grey60", size=0.2) +
geom_line( aes(x=YearQuarter, y=Mean), col = "red") +
facet_wrap(vars(Population), ncol = 1) +
scale_x_date(date_labels = "%Y", date_breaks = "1 year") +
labs(x = "Year", y = "Number of services per 1000 population") +
theme_bw() +
theme(
strip.background=element_rect(fill="white")
)
IRSAD distribution
IRSAD<-dta_imputed %>%
dplyr::select(SA3_NAME,IRSAD) %>%
distinct()
summary(IRSAD)## SA3_NAME IRSAD
## Adelaide City : 1 Min. : 748.6
## Adelaide Hills: 1 1st Qu.: 943.3
## Albany : 1 Median : 983.1
## Albury : 1 Mean : 995.5
## Alice Springs : 1 3rd Qu.:1044.8
## Armadale : 1 Max. :1166.6
## (Other) :325a<-IRSAD %>%
ggplot(aes(IRSAD)) +
stat_ecdf(geom = "step") +
theme_bw() +
scale_y_continuous(labels=scales::percent)+
scale_x_continuous(limits = c(740,1170)) +
labs(y="Percentage")
b<-IRSAD %>%
mutate(IRSAD=plyr::round_any(IRSAD,30)) %>%
ggplot(aes(IRSAD)) +
geom_bar(aes(y = (..count..)/sum(..count..))) +
theme_bw() +
scale_y_continuous(labels=scales::percent) +
scale_x_continuous(limits = c(740,1170)) +
labs(y="Percentage")
ggpubr::ggarrange(a, b,
labels = c("A", "B"),nrow = 2)
# Descriptive table
This table describes the median (IQR) number of services delivered per 1000 people at an SA3 level for 2019 and 2020 for each subgroup, as well as the percentage change in these values.
tbl_dta<-dta %>%
filter(Year>=2019) %>%
group_by(SA3_CODE16,Population,Year) %>%
summarise(Services=sum(Services), ERP=mean(ERP)) %>%
mutate(Rates=Services*1000/ERP) %>%
ungroup() %>%
dplyr::select(Population,Year,Rates,SA3_CODE16) %>%
pivot_wider(names_from=Year,values_from=Rates) %>%
mutate(`Percentage increase`= ( `2020`- `2019`)*100/ `2019`)
tbl_dta %>%
dplyr::select(-SA3_CODE16) %>%
tbl_summary(
by = Population,
missing = "no",
digits = `Percentage increase` ~ 1,
) | Characteristic | Age 12-17, N = 3321 | Age 18-25, N = 3321 | Females, N = 3321 | Males, N = 3321 | All young people, N = 3321 |
|---|---|---|---|---|---|
| 2019 | 572 (428, 688) | 696 (552, 850) | 837 (661, 1,008) | 438 (346, 543) | 639 (496, 760) |
| 2020 | 622 (470, 787) | 793 (618, 997) | 984 (779, 1,231) | 443 (360, 569) | 717 (551, 889) |
| Percentage increase | 12.6 (4.0, 21.0) | 15.7 (9.9, 21.8) | 20.4 (14.0, 28.0) | 2.8 (-2.9, 8.5) | 14.0 (8.4, 20.3) |
|
1
Median (IQR)
|
|||||
Boxplot
These boxplots show the distribution of the number of services accessed per 1000 people for each quarter between 2019 and 2020 by population subgroup and year.
dta %>%
filter(Year>=2019) %>%
mutate(Quarter=as.factor(paste0("Q", Quarter)),
Year=as.factor(Year)) %>%
mutate(Rates=Services*1000/ERP) %>%
ggplot(aes(x=Quarter, y=Rates, col=Year)) +
geom_boxplot() +theme_bw() +
facet_wrap( ~ Population,ncol = 5)+
theme(legend.position = "bottom",
strip.background=element_rect(fill="white"))+
scale_color_manual(values=c("#0072B2","#E69F00")) +
labs(y="Rates per 1000 population")
dta %>%
filter(Year>=2019) %>%
mutate(Quarter=as.factor(paste0("Q", Quarter)),
Year=as.factor(Year)) %>%
mutate(Rates=Services*1000/ERP) %>%
dplyr::select(Quarter,Rates,Year,Population,SA3_NAME) %>%
pivot_wider(names_from="Population",values_from = "Rates") %>%
dplyr::select(-SA3_NAME) %>%
tbl_strata(
strata = Quarter,
.tbl_fun =
~ .x %>%
tbl_summary(by = Year, missing = "no")
)| Characteristic | Q1 | Q2 | Q3 | Q4 | ||||
|---|---|---|---|---|---|---|---|---|
| 2019, N = 3321 | 2020, N = 3321 | 2019, N = 3321 | 2020, N = 3321 | 2019, N = 3321 | 2020, N = 3321 | 2019, N = 3321 | 2020, N = 3321 | |
| Age 12-17 | 132 (98, 161) | 134 (100, 168) | 148 (111, 181) | 148 (111, 192) | 152 (114, 185) | 173 (132, 226) | 134 (99, 165) | 163 (119, 211) |
| Age 18-25 | 167 (133, 206) | 176 (137, 219) | 173 (136, 215) | 197 (152, 245) | 183 (142, 223) | 219 (169, 272) | 169 (131, 204) | 202 (158, 250) |
| Females | 199 (155, 244) | 209 (163, 260) | 214 (168, 260) | 243 (187, 306) | 223 (173, 264) | 281 (216, 343) | 198 (155, 239) | 260 (203, 321) |
| Males | 108 (81, 132) | 106 (82, 131) | 111 (88, 138) | 110 (90, 140) | 118 (91, 145) | 122 (97, 156) | 109 (85, 131) | 108 (85, 140) |
| All young people | 151 (117, 184) | 155 (122, 195) | 161 (126, 198) | 174 (135, 219) | 168 (131, 203) | 200 (155, 249) | 153 (119, 184) | 180 (142, 230) |
|
1
Median (IQR)
|
||||||||
State comparision
Conversely, these boxplots depict the distribution of the yearly number of services delivered per 1000 people for Victoria and the rest of Australia by subpopulation and year. The yearly rates are obtained similarly to the formula for the heat maps above.
dta %>%
filter(Year>=2019 ) %>%
mutate(Quarter=as.factor(paste0("Q", Quarter)),
Year=as.factor(Year)) %>%
group_by(SA3_CODE16,Population,Year,State) %>%
summarise(Services=sum(Services), ERP=mean(ERP)) %>%
mutate(Rates=Services*1000/ERP) %>%
ungroup() %>%
ggplot(aes(x=State, y=Rates, col=Year)) +
geom_boxplot() +theme_bw() +
facet_wrap( ~ Population,ncol = 5)+
theme(legend.position = "bottom",
strip.background=element_rect(fill="white"))+
scale_color_manual(values=c("#0072B2","#E69F00")) +
labs(y="Rates per 1000 population")
dta %>%
filter(Year>=2019 ) %>%
mutate(Quarter=as.factor(paste0("Q", Quarter)),
Year=as.factor(Year)) %>%
group_by(SA3_CODE16,Population,Year,State) %>%
summarise(Services=sum(Services), ERP=mean(ERP)) %>%
mutate(Rates=Services*1000/ERP) %>%
ungroup() %>%
dplyr::select(State,Rates,Year,Population,SA3_CODE16) %>%
pivot_wider(names_from="Year",values_from = "Rates") %>%
mutate(`Percentage increase`= (`2020`-`2019`)*100/`2019`) %>%
dplyr::select(-SA3_CODE16) %>%
tbl_strata(
strata = Population,
.tbl_fun =
~ .x %>%
tbl_summary(by = State, missing = "no",
digits = `Percentage increase` ~ 1)
)| Characteristic | Age 12-17 | Age 18-25 | Females | Males | All young people | |||||
|---|---|---|---|---|---|---|---|---|---|---|
| Other, N = 2661 | Victoria, N = 661 | Other, N = 2661 | Victoria, N = 661 | Other, N = 2661 | Victoria, N = 661 | Other, N = 2661 | Victoria, N = 661 | Other, N = 2661 | Victoria, N = 661 | |
| 2019 | 553 (413, 680) | 636 (508, 733) | 682 (525, 839) | 765 (610, 870) | 813 (639, 999) | 921 (751, 1,089) | 435 (335, 536) | 474 (405, 550) | 612 (482, 757) | 695 (560, 817) |
| 2020 | 619 (457, 783) | 655 (511, 855) | 761 (595, 972) | 886 (739, 1,071) | 963 (744, 1,205) | 1,155 (933, 1,405) | 439 (344, 556) | 465 (408, 596) | 685 (540, 874) | 796 (659, 984) |
| Percentage increase | 13.3 (5.3, 21.3) | 8.9 (2.1, 17.4) | 14.3 (8.8, 20.2) | 20.9 (15.3, 26.2) | 19.4 (13.4, 27.1) | 25.6 (16.3, 31.1) | 2.9 (-3.2, 8.6) | 2.5 (-2.7, 7.7) | 13.1 (8.2, 19.8) | 16.9 (10.1, 22.3) |
|
1
Median (IQR)
|
||||||||||
Using data from 2015-2020
Indiviudal model
The interrupted time series models for the mean number of services at time \(t\), \(\mu_t\), took the form:
\[\log(\mu_t) = \beta_0 + \beta_1t + \beta_2COVID_t + \beta_{3}Quarter_{t} + \log(ERP_t)\]
,where \(COVID_t\) is the interruption jump coefficient, \(Quarter_{t}\) is the quarterly seasonal effect term and \(ERP_t\) is the estimated resident population at time \(t\).
#initialize results
results<-list()
count=1
# loop indiviudal time sereis
for(id in unique(sort(dta_imputed$ID))) {
temp <-dta_imputed %>%
filter(ID==id & Year >= 2015)
#check for evidence of residual autocorrelation
model <- glm(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)) ,
data = temp, family = "quasipoisson")
res<-resid(model, type='pearson')
test<-Box.test(res, lag=10, fitdf=0)
# when p is large will not control for autocorrelation
if(test$p.value>0.05){
#convergence issues with the default optimiser in a few series
if(tryCatch(model<-MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)),
random = ~ 1|SA3_NAME,
data =temp, family = "quasipoisson"),
error = function(c) "error")=="error")
{
model <- MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)),
random = ~ 1|SA3_NAME, temp, family = "quasipoisson",
control = lmeControl(opt = "optim"))
}
} else {
count=count+1
#convergence issues with the default optimiser in a few series
if(tryCatch(model <- MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)),
random = ~ 1|SA3_NAME, data = temp, family = "quasipoisson",
correlation = corAR1()), error = function(c) "error")=="error") {
model <- MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)),
random = ~ 1|SA3_NAME, data = temp, family = "quasipoisson",
control = lmeControl(opt = "optim"),
correlation = corAR1())
}
}
#save results
result <- coef(summary(model))
result <-cbind(rownames(result),as_tibble(result))
result$SA3_CODE16 <-unique(temp[temp$ID==id,]$SA3_CODE16)
result$Population <-unique(temp[temp$ID==id,]$Population)
results[[id]]<-result
}
#extract SA3 information
sa3 <- dta %>%
dplyr::select(SA3_CODE16,SA3_NAME,IRSAD,Region, State) %>%
unique()
results<-as_tibble(rbindlist(results,idcol="ID")) %>%
mutate(SA3_CODE16=as.numeric(SA3_CODE16)) %>%
left_join(sa3,by="SA3_CODE16")
names(results)[2]<-"Terms"
#fraction of models with autocorrelation issues
count/length(unique(sort(dta_imputed$ID)))## [1] 0.1212675write.csv(results %>% filter(Terms=="COVID"),"results.csv")Plots
This scatter plot shows the estimated relative risk against the IRSAD score for each subgroup. Superimposed is locally weighted fitted smoothed line and the 95% CI for the prediction at that IRSAD score.
COVID<-results %>%
filter(Terms=="COVID") %>%
mutate(RR=exp(Value))
COVID %>%
#3 sa3 excluded at the tail
dplyr::filter( IRSAD > 810) %>%
ggplot(aes(x=IRSAD, y=RR)) +
geom_point( alpha=0.6, col="#0072B2", size=1)+
geom_smooth(col="red",span=1) +
facet_wrap( ~ Population,ncol = 5)+
theme_bw() +
facet_wrap( ~ Population,ncol = 5)+
theme(legend.position = "bottom",
strip.background=element_rect(fill="white"))
These boxplots compare the distribution of the estimated relative risks by state (Victoria vs. rest of Australia) and region (metropolitan vs. reginal) for each subpopulation.
state<-COVID %>%
ggplot(aes(x=State, y=RR, colour=State)) +
geom_boxplot()+
facet_wrap( ~ Population,ncol = 5)+
theme_bw() +
facet_wrap( ~ Population,ncol = 5)+
theme(legend.position = "none",
strip.background=element_rect(fill="white"))+
scale_color_manual(values=c("#56B4E9","#0072B2"))
region<- COVID %>%
ggplot(aes(x=Region, y=RR, colour=Region)) +
geom_boxplot()+
facet_wrap( ~ Population,ncol = 5)+
theme_bw() +
facet_wrap( ~ Population,ncol = 5)+
theme(legend.position = "none",
strip.background=element_rect(fill="white"))+
scale_color_manual(values=c("#E69F00", "#D55E00"))
ggpubr::ggarrange(state,region,
labels = c("A", "B"), nrow = 2)
These two scatter plots compare the estimated relative risks within each population variable (age group and sex) and their relationships with the colouring variable, IRSAD score. The dashed line represents equality between the estimates for each of the sexes/age-groups.
library(ggpubr)
library(wesanderson)
pal <- wes_palette("Zissou1", 100, type = "continuous")
a<-COVID %>%
dplyr::select(SA3_CODE16,Population,RR,IRSAD) %>%
filter(Population=="Males" | Population=="Females") %>%
pivot_wider(names_from=Population,values_from=RR) %>%
mutate(compare=Females>Males) %>%
ggplot(aes(x=Males, y=Females, colour=IRSAD)) +
geom_jitter( alpha=0.8, size=2)+
theme_bw() +
scale_color_gradientn(colours = pal) +
scale_x_continuous(limits = c(0.6, 1.4)) +
scale_y_continuous(limits = c(0.6, 1.4)) +
geom_abline(slope = 1, col="gray50",linetype = 2)+
labs(x="RR for males", y="RR for females") +
theme(legend.position="none")
b<-COVID %>%
dplyr::select(SA3_CODE16,Population,RR,IRSAD) %>%
filter(Population=="Age 12-17" | Population=="Age 18-25") %>%
pivot_wider(names_from=Population,values_from=RR) %>%
mutate(compare=`Age 18-25` >`Age 12-17`) %>%
ggplot(aes(x=`Age 18-25`, y=`Age 12-17`, colour=IRSAD)) +
geom_jitter( alpha=0.8, size=2)+
theme_bw() +
scale_color_gradientn(colours = pal) +
scale_x_continuous(limits = c(0.6, 1.4)) +
scale_y_continuous(limits = c(0.6, 1.4)) +
geom_abline(slope = 1, col="gray50",linetype = 2)+
labs(x="RR for age 12-17", y="RR for age 18-25")+
theme(legend.position="none")
g_legend<-function(a.gplot){
tmp <- ggplot_gtable(ggplot_build(a.gplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
legend<-a +
theme(legend.position="bottom") +
theme(legend.key.height= unit(0.5, 'cm'),
legend.key.width= unit(1.2, 'cm'))
legend <- g_legend(legend)
grid.arrange(ggarrange(a,b,ncol=2) ,
legend, nrow=2,heights=c(10, 1))
Pooled effect
This forest plot portrays the average estimated relative risk (95% CI) for each subgroup as estimated by a random effects meta-analysis model.
library(metafor)
extract_meta<-function(data){
a<-rma(Value, Std.Error, data=data,method = "REML", test = "knha")
ret<-exp(c(a$b,a$ci.lb, a$ci.ub ))
names(ret)<-c( "RR", "CIL","CIH")
ret
}
res <- results %>%
filter(Terms=="COVID") %>%
group_by(Population) %>%
group_modify(~ broom::tidy(extract_meta(data = .x))) %>%
pivot_wider(names_from=names,values_from=x) %>%
mutate('RR (95% CI)' = paste0(format(round(RR, digits=3), nsmall=3), " (" ,
format(round( CIL, digits = 3),nsmall=3)," - " ,
format(round(CIH, digits = 3),nsmall=3),")"))
res$y_loc=5:1
#main plot in the middle
plot1<- ggplot(res, aes(x = RR, y = y_loc, colour ) )+
geom_point() +
geom_errorbarh(aes(xmin = CIL, xmax = CIH), height = 0.3) +
ylab(NULL) +
scale_y_continuous(limits = c(0.8,5.5)) +
ggtitle(" ") +
xlab("Estimated RR (95% CI)") +
scale_colour_manual(values = colours ) +
theme(plot.title = element_text(hjust =1 ,size = 12,face = "bold"),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.border = element_blank(),
axis.line.x = element_line(colour = "black"),
panel.background = element_blank(),
panel.grid = element_blank(),
plot.margin = margin(0, 0, 0, 0, "cm"),
legend.position = "none") +
geom_vline(aes(xintercept = 1), linetype = "dashed", col="gray50")
#base plot for adding text
tab_base <- ggplot(res, aes(y = y_loc)) +
ylab(NULL) + xlab("") +
scale_y_continuous(limits = c(0.8,5.5)) +
theme(plot.title = element_text( size = 12), ## centering title on text
axis.text.x = element_text(color = "white"), ## need text to be printed so it stays aligned with figure but white so it's invisible
axis.line = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.y = element_blank(),
plot.margin = margin(0, 0, 0.3, 0, "cm"),
legend.position = "none",
panel.background = element_blank(),
panel.grid = element_blank()) +
scale_colour_manual(values = colours ) +
labs(x=" ")
# Predictors
tab1 <- tab_base +
geom_text(aes(x = 0, label = Population, y = y_loc), hjust = 0,nudge_x = 0.05) +
scale_x_continuous(limits = c(0,0.4)) +
ggtitle("Outcomes") +
theme(plot.title = element_text(hjust = 0.3))
# OR
OR1 <- tab_base +
geom_text(aes(x = 0, label = `RR (95% CI)`) , hjust = 0.5) +
ggtitle(" OR (95% CI)")+
theme(plot.title = element_text(hjust = 0.5))
gridExtra::grid.arrange(tab1, plot1, OR1,
layout_matrix = matrix(c(rep(1,2),
rep(2,6),
rep(3,3)), nrow = 1))
Meta regression
Random effects meta-regression was used to compile the interruption effect estimates and identify their relationships with the predictors IRSAD score, state (rest of Australia vs. Victoria) and area type (metropolitan vs. regional). The meta-regression model took the form:
\[\hat{\beta}^{i}_2 = \alpha_0 + \alpha_1Regional_i+\alpha_2IRSAD_i+ \alpha_3Victoria_i + u_{i} + \epsilon_{i}\]
, where \(\hat{\beta}^{ij}_2\) is the estimated coefficient for the COVID-19 impact for SA3 \(i\) in a specific population. \(Regional_i\) is the binary area type variable (baseline metropolitan), \(IRSAD_i\) is the IRSAD score for SA3 \(i\) centred on the mean IRSAD score for Australia and expressed in 100 units, \(Victoria_i\) is the binary state variable (baseline rest of Australia), \(\epsilon_{i}\) is the sampling error and \(u_{i}\) is the error due to heterogeneity between SA3s in true interruption effect size.
library(metafor)
reg_meta<-function(data){
a<-rma(yi=Value, sei=Std.Error, data=data,method = "REML",
mods = ~ IRSAD + State + Region, test = "knha")
browser()
print( unique(data$population))
print(summary(a))
ret<-data.frame(cbind(exp(cbind(a$b, a$ci.lb,a$ci.ub)), a$pval))
names(ret)<-c( "RR", "CIL","CIH", "p")
ret$Variable<-c("Intercept","IRSAD","Victoria", "Regional" )
return(ret)
}
res <- results %>%
filter(Terms=="COVID") %>%
mutate(IRSAD= (IRSAD-995.5)/100) %>%
mutate(population=Population)%>%
group_by(Population) %>%
group_modify(~reg_meta(data = .x)) %>%
mutate('RR ' = format(round(RR, digits=3), nsmall=3),
"95%CI" =paste0(format(round( CIL, digits = 3),nsmall=3)," - " ,
format(round(CIH, digits = 3),nsmall=3)),
'p-value'=ifelse(p<=0.001, "<0.001", format(round(p, digits=3), nsmall=3))) %>%
dplyr::select(-RR,-CIL,-CIH,-p) %>%
pivot_longer(names_to="Estimate", values_to="Value",-c("Population","Variable") ) %>%
pivot_wider(names_from = Population,values_from=Value)## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## [1] Age 12-17
## Levels: Age 12-17 Age 18-25 Females Males All young people
## debug at <text>#7: print(summary(a))
##
## Mixed-Effects Model (k = 327; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 250.9328 -501.8655 -491.8655 -472.9773 -491.6762
##
## tau^2 (estimated amount of residual heterogeneity): 0.0071 (SE = 0.0009)
## tau (square root of estimated tau^2 value): 0.0840
## I^2 (residual heterogeneity / unaccounted variability): 67.04%
## H^2 (unaccounted variability / sampling variability): 3.03
## R^2 (amount of heterogeneity accounted for): 20.88%
##
## Test for Residual Heterogeneity:
## QE(df = 323) = 971.5387, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 323) = 19.5891, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.0337 0.0090 3.7440 0.0002 0.0160 0.0515 ***
## IRSAD 0.0746 0.0105 7.1266 <.0001 0.0540 0.0952 ***
## StateVictoria -0.0240 0.0145 -1.6615 0.0976 -0.0525 0.0044 .
## RegionRegional 0.0223 0.0143 1.5534 0.1213 -0.0059 0.0505
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## [1] Age 18-25
## Levels: Age 12-17 Age 18-25 Females Males All young people
## debug at <text>#7: print(summary(a))
##
## Mixed-Effects Model (k = 328; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 354.9734 -709.9469 -699.9469 -681.0432 -699.7582
##
## tau^2 (estimated amount of residual heterogeneity): 0.0030 (SE = 0.0004)
## tau (square root of estimated tau^2 value): 0.0550
## I^2 (residual heterogeneity / unaccounted variability): 59.08%
## H^2 (unaccounted variability / sampling variability): 2.44
## R^2 (amount of heterogeneity accounted for): 44.50%
##
## Test for Residual Heterogeneity:
## QE(df = 324) = 800.0453, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 324) = 41.8378, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.0374 0.0064 5.8261 <.0001 0.0248 0.0500 ***
## IRSAD 0.0495 0.0074 6.6672 <.0001 0.0349 0.0641 ***
## StateVictoria 0.0962 0.0109 8.8435 <.0001 0.0748 0.1176 ***
## RegionRegional 0.0304 0.0103 2.9648 0.0033 0.0102 0.0506 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## [1] Females
## Levels: Age 12-17 Age 18-25 Females Males All young people
## debug at <text>#7: print(summary(a))
##
## Mixed-Effects Model (k = 327; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 346.9755 -693.9511 -683.9511 -665.0628 -683.7618
##
## tau^2 (estimated amount of residual heterogeneity): 0.0029 (SE = 0.0004)
## tau (square root of estimated tau^2 value): 0.0542
## I^2 (residual heterogeneity / unaccounted variability): 54.42%
## H^2 (unaccounted variability / sampling variability): 2.19
## R^2 (amount of heterogeneity accounted for): 42.61%
##
## Test for Residual Heterogeneity:
## QE(df = 323) = 716.2943, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 323) = 37.7070, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.0833 0.0065 12.7335 <.0001 0.0705 0.0962 ***
## IRSAD 0.0664 0.0077 8.6598 <.0001 0.0513 0.0815 ***
## StateVictoria 0.0636 0.0111 5.7289 <.0001 0.0418 0.0854 ***
## RegionRegional 0.0339 0.0105 3.2279 0.0014 0.0132 0.0546 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## [1] Males
## Levels: Age 12-17 Age 18-25 Females Males All young people
## debug at <text>#7: print(summary(a))
##
## Mixed-Effects Model (k = 328; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 293.9581 -587.9161 -577.9161 -559.0124 -577.7274
##
## tau^2 (estimated amount of residual heterogeneity): 0.0044 (SE = 0.0006)
## tau (square root of estimated tau^2 value): 0.0667
## I^2 (residual heterogeneity / unaccounted variability): 63.89%
## H^2 (unaccounted variability / sampling variability): 2.77
## R^2 (amount of heterogeneity accounted for): 11.88%
##
## Test for Residual Heterogeneity:
## QE(df = 324) = 892.8698, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 324) = 9.7570, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt -0.0494 0.0078 -6.3100 <.0001 -0.0648 -0.0340 ***
## IRSAD 0.0449 0.0089 5.0314 <.0001 0.0273 0.0624 ***
## StateVictoria 0.0141 0.0125 1.1269 0.2606 -0.0105 0.0386
## RegionRegional 0.0165 0.0124 1.3292 0.1847 -0.0079 0.0410
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## [1] All young people
## Levels: Age 12-17 Age 18-25 Females Males All young people
## debug at <text>#7: print(summary(a))
##
## Mixed-Effects Model (k = 331; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 382.8617 -765.7234 -755.7234 -736.7736 -755.5364
##
## tau^2 (estimated amount of residual heterogeneity): 0.0026 (SE = 0.0004)
## tau (square root of estimated tau^2 value): 0.0515
## I^2 (residual heterogeneity / unaccounted variability): 60.75%
## H^2 (unaccounted variability / sampling variability): 2.55
## R^2 (amount of heterogeneity accounted for): 39.45%
##
## Test for Residual Heterogeneity:
## QE(df = 327) = 835.6417, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 327) = 35.2839, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.0373 0.0059 6.2989 <.0001 0.0256 0.0489 ***
## IRSAD 0.0591 0.0069 8.6269 <.0001 0.0457 0.0726 ***
## StateVictoria 0.0498 0.0100 4.9931 <.0001 0.0302 0.0694 ***
## RegionRegional 0.0266 0.0095 2.8057 0.0053 0.0079 0.0452 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)res$Variable[c(2:3,5:6,8:9,11:12)]<-""
print_table(res, align = c("l", "l",rep("c",5)), colwidth = c(2,3,rep(4.4,5)))## <caption> (\#tab:)</caption>Variable | Estimate | Age 12-17 | Age 18-25 | Females | Males | All young people |
Intercept | RR | 1.034 | 1.038 | 1.087 | 0.952 | 1.038 |
95%CI | 1.016 - 1.053 | 1.025 - 1.051 | 1.073 - 1.101 | 0.937 - 0.967 | 1.026 - 1.050 | |
p-value | <0.001 | <0.001 | <0.001 | <0.001 | <0.001 | |
IRSAD | RR | 1.077 | 1.051 | 1.069 | 1.046 | 1.061 |
95%CI | 1.055 - 1.100 | 1.035 - 1.066 | 1.053 - 1.085 | 1.028 - 1.064 | 1.047 - 1.075 | |
p-value | <0.001 | <0.001 | <0.001 | <0.001 | <0.001 | |
Victoria | RR | 0.976 | 1.101 | 1.066 | 1.014 | 1.051 |
95%CI | 0.949 - 1.004 | 1.078 - 1.125 | 1.043 - 1.089 | 0.990 - 1.039 | 1.031 - 1.072 | |
p-value | 0.098 | <0.001 | <0.001 | 0.261 | <0.001 | |
Regional | RR | 1.023 | 1.031 | 1.034 | 1.017 | 1.027 |
95%CI | 0.994 - 1.052 | 1.010 - 1.052 | 1.013 - 1.056 | 0.992 - 1.042 | 1.008 - 1.046 | |
p-value | 0.121 | 0.003 | 0.001 | 0.185 | 0.005 |
Sensitivity analysis using data from 2011
All results below are generated similarly as before for longer time series from 2011.
Individual model
#initialize results
results<-list()
count=1
# loop indiviudal time sereis
for(id in unique(sort(dta_imputed$ID))) {
temp <-dta_imputed %>%
filter(ID==id )
#check for evidence of residual autocorrelation
model <- glm(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)) ,
data = temp, family = "quasipoisson")
res<-resid(model, type='pearson')
test<-Box.test(res, lag=10, fitdf=0)
# when p is large will not control for autocorrelation
if(test$p.value>0.05){
#convergence issues with the default optimiser in a few series
if(tryCatch(model<-MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)),
random = ~ 1|SA3_NAME,
data =temp, family = "quasipoisson"),
error = function(c) "error")=="error")
{
model <- MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)),
random = ~ 1|SA3_NAME, temp, family = "quasipoisson",
control = lmeControl(opt = "optim"))
}
} else {
count=count+1
#convergence issues with the default optimiser in a few series
if(tryCatch(model <- MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)),
random = ~ 1|SA3_NAME, data = temp, family = "quasipoisson",
correlation = corAR1()), error = function(c) "error")=="error") {
model <- MASS::glmmPQL(Services ~ YearQuarter + COVID + Quarter + offset(log(ERP)),
random = ~ 1|SA3_NAME, data = temp, family = "quasipoisson",
control = lmeControl(opt = "optim"),
correlation = corAR1())
}
}
#save results
result <- coef(summary(model))
result <-cbind(rownames(result),as_tibble(result))
result$SA3_CODE16 <-unique(temp[temp$ID==id,]$SA3_CODE16)
result$Population <-unique(temp[temp$ID==id,]$Population)
results[[id]]<-result
}
#extract SA3 information
sa3 <- dta %>%
dplyr::select(SA3_CODE16,SA3_NAME,IRSAD,Region, State) %>%
unique()
results<-as_tibble(rbindlist(results,idcol="ID")) %>%
mutate(SA3_CODE16=as.numeric(SA3_CODE16)) %>%
left_join(sa3,by="SA3_CODE16")
names(results)[2]<-"Terms"
#fraction of models with autocorrelation issues
count/length(unique(sort(dta_imputed$ID)))## [1] 0.4582572Plots
COVID<-results %>%
filter(Terms=="COVID") %>%
mutate(RR=exp(Value))
COVID %>%
#3 sa3 excluded at the tail
dplyr::filter( IRSAD > 810 ) %>%
ggplot(aes(x=IRSAD, y=RR)) +
geom_point( alpha=0.6, col="#0072B2", size=1)+
geom_smooth(col="red",span=1) +
facet_wrap( ~ Population,ncol = 5)+
theme_bw() +
facet_wrap( ~ Population,ncol = 5)+
theme(legend.position = "bottom",
strip.background=element_rect(fill="white"))
state<-COVID %>%
ggplot(aes(x=State, y=RR, colour=State)) +
geom_boxplot()+
facet_wrap( ~ Population,ncol = 5)+
theme_bw() +
facet_wrap( ~ Population,ncol = 5)+
theme(legend.position = "none",
strip.background=element_rect(fill="white"))+
scale_color_manual(values=c("#56B4E9","#0072B2"))
region<- COVID %>%
ggplot(aes(x=Region, y=RR, colour=Region)) +
geom_boxplot()+
facet_wrap( ~ Population,ncol = 5)+
theme_bw() +
facet_wrap( ~ Population,ncol = 5)+
theme(legend.position = "none",
strip.background=element_rect(fill="white"))+
scale_color_manual(values=c("#E69F00", "#D55E00"))
ggpubr::ggarrange(state,region,
labels = c("A", "B"), nrow = 2)
library(ggpubr)
library(wesanderson)
pal <- wes_palette("Zissou1", 100, type = "continuous")
a<-COVID %>%
dplyr::select(SA3_CODE16,Population,RR,IRSAD) %>%
filter(Population=="Males" | Population=="Females") %>%
pivot_wider(names_from=Population,values_from=RR) %>%
mutate(compare=Females>Males) %>%
ggplot(aes(x=Males, y=Females, colour=IRSAD)) +
geom_jitter( alpha=0.8, size=2)+
theme_bw() +
scale_color_gradientn(colours = pal) +
scale_x_continuous(limits = c(0.6, 1.4)) +
scale_y_continuous(limits = c(0.6, 1.4)) +
geom_abline(slope = 1, col="gray50",linetype = 2)+
labs(x="RR for males", y="RR for females") +
theme(legend.position="none")
b<-COVID %>%
dplyr::select(SA3_CODE16,Population,RR,IRSAD) %>%
filter(Population=="Age 12-17" | Population=="Age 18-25") %>%
pivot_wider(names_from=Population,values_from=RR) %>%
mutate(compare=`Age 18-25` >`Age 12-17`) %>%
ggplot(aes(x=`Age 18-25`, y=`Age 12-17`, colour=IRSAD)) +
geom_jitter( alpha=0.8, size=2)+
theme_bw() +
scale_color_gradientn(colours = pal) +
scale_x_continuous(limits = c(0.6, 1.4)) +
scale_y_continuous(limits = c(0.6, 1.4)) +
geom_abline(slope = 1, col="gray50",linetype = 2)+
labs(x="RR for age 12-17", y="RR for age 18-25")+
theme(legend.position="none")
g_legend<-function(a.gplot){
tmp <- ggplot_gtable(ggplot_build(a.gplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
legend<-a +
theme(legend.position="bottom") +
theme(legend.key.height= unit(0.5, 'cm'),
legend.key.width= unit(1.2, 'cm'))
legend <- g_legend(legend)
grid.arrange(ggarrange(a,b,ncol=2) ,
legend, nrow=2,heights=c(10, 1))
Pooled effect
res <- results %>%
filter(Terms=="COVID") %>%
group_by(Population) %>%
group_modify(~ broom::tidy(extract_meta(data = .x))) %>%
pivot_wider(names_from=names,values_from=x) %>%
mutate('RR (95% CI)' = paste0(format(round(RR, digits=3), nsmall=3), " (" ,
format(round( CIL, digits = 3),nsmall=3)," - " ,
format(round(CIH, digits = 3),nsmall=3),")"))
res$y_loc=5:1
#main plot in the middle
plot1<- ggplot(res, aes(x = RR, y = y_loc, colour ) )+
geom_point() +
geom_errorbarh(aes(xmin = CIL, xmax = CIH), height = 0.3) +
ylab(NULL) +
scale_y_continuous(limits = c(0.8,5.5)) +
ggtitle(" ") +
xlab("Estimated RR (95% CI)") +
scale_colour_manual(values = colours ) +
theme(plot.title = element_text(hjust =1 ,size = 12,face = "bold"),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.border = element_blank(),
axis.line.x = element_line(colour = "black"),
panel.background = element_blank(),
panel.grid = element_blank(),
plot.margin = margin(0, 0, 0, 0, "cm"),
legend.position = "none") +
geom_vline(aes(xintercept = 1), linetype = "dashed", col="gray50")
#base plot for adding text
tab_base <- ggplot(res, aes(y = y_loc)) +
ylab(NULL) + xlab("") +
scale_y_continuous(limits = c(0.8,5.5)) +
theme(plot.title = element_text( size = 12), ## centering title on text
axis.text.x = element_text(color = "white"), ## need text to be printed so it stays aligned with figure but white so it's invisible
axis.line = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.y = element_blank(),
plot.margin = margin(0, 0, 0.3, 0, "cm"),
legend.position = "none",
panel.background = element_blank(),
panel.grid = element_blank()) +
scale_colour_manual(values = colours ) +
labs(x=" ")
# Predictors
tab1 <- tab_base +
geom_text(aes(x = 0, label = Population, y = y_loc), hjust = 0,nudge_x = 0.05) +
scale_x_continuous(limits = c(0,0.4)) +
ggtitle("Outcomes") +
theme(plot.title = element_text(hjust = 0.3))
# OR
OR1 <- tab_base +
geom_text(aes(x = 0, label = `RR (95% CI)`) , hjust = 0.5) +
ggtitle(" OR (95% CI)")+
theme(plot.title = element_text(hjust = 0.5))
gridExtra::grid.arrange(tab1, plot1, OR1,
layout_matrix = matrix(c(rep(1,2),
rep(2,6),
rep(3,3)), nrow = 1))
Meta regression
res <- results %>%
filter(Terms=="COVID") %>%
mutate(IRSAD= (IRSAD-995.5)/100) %>%
group_by(Population) %>%
group_modify(~reg_meta(data = .x)) %>%
mutate('RR ' = format(round(RR, digits=3), nsmall=3),
"95%CI" =paste0(format(round( CIL, digits = 3),nsmall=3)," - " ,
format(round(CIH, digits = 3),nsmall=3)),
'p-value'=ifelse(p<=0.001, "<0.001", format(round(p, digits=3), nsmall=3))) %>%
dplyr::select(-RR,-CIL,-CIH,-p) %>%
pivot_longer(names_to="Estimate", values_to="Value",-c("Population","Variable") ) %>%
pivot_wider(names_from = Population,values_from=Value)
res$Variable[c(2:3,5:6,8:9,11:12)]<-""
print_table(res, align = c("l", "l",rep("c",5)), colwidth = c(2,3,rep(4.4,5)))
Sensitivity analysis exclude autocorrelation
All results below are generated similarly as before for interrupted time series model without AR(1) structure estimated using the glm function.
Individual model
results<-list()
for(id in unique(sort(dta_imputed$ID))) {
temp <-dta_imputed %>%
filter(ID==id & Year >= 2015)
model <- glm(Services ~ Year + COVID + Quarter + offset(log(ERP)) , data = temp, family = "quasipoisson")
result <- coef(summary(model))
result <-cbind(rownames(result),as_tibble(result))
result$SA3_CODE16 <-unique(temp[temp$ID==id,]$SA3_CODE16)
result$Population <-unique(temp[temp$ID==id,]$Population)
results[[id]]<-result
}
#combine results
results<-as_tibble(rbindlist(results,idcol="ID")) %>%
mutate(SA3_CODE16=as.numeric(SA3_CODE16)) %>%
left_join(sa3,by="SA3_CODE16")
names(results)[2:4]<-c("Terms", "Value","Std.Error")Plots
COVID<-results %>%
filter(Terms=="COVID") %>%
mutate(RR=exp(Value))
COVID %>%
#3 sa3 excluded at the tail
dplyr::filter( IRSAD > 810 ) %>%
ggplot(aes(x=IRSAD, y=RR)) +
geom_point( alpha=0.6, col="#0072B2", size=1)+
geom_smooth(col="red",span=1) +
facet_wrap( ~ Population,ncol = 5)+
theme_bw() +
facet_wrap( ~ Population,ncol = 5)+
theme(legend.position = "bottom",
strip.background=element_rect(fill="white"))
state<-COVID %>%
ggplot(aes(x=State, y=RR, colour=State)) +
geom_boxplot()+
facet_wrap( ~ Population,ncol = 5)+
theme_bw() +
facet_wrap( ~ Population,ncol = 5)+
theme(legend.position = "none",
strip.background=element_rect(fill="white"))+
scale_color_manual(values=c("#56B4E9","#0072B2"))
region<- COVID %>%
ggplot(aes(x=Region, y=RR, colour=Region)) +
geom_boxplot()+
facet_wrap( ~ Population,ncol = 5)+
theme_bw() +
facet_wrap( ~ Population,ncol = 5)+
theme(legend.position = "none",
strip.background=element_rect(fill="white"))+
scale_color_manual(values=c("#E69F00", "#D55E00"))
ggpubr::ggarrange(state,region,
labels = c("A", "B"), nrow = 2)
library(ggpubr)
library(wesanderson)
pal <- wes_palette("Zissou1", 100, type = "continuous")
a<-COVID %>%
dplyr::select(SA3_CODE16,Population,RR,IRSAD) %>%
filter(Population=="Males" | Population=="Females") %>%
pivot_wider(names_from=Population,values_from=RR) %>%
mutate(compare=Females>Males) %>%
ggplot(aes(x=Males, y=Females, colour=IRSAD)) +
geom_jitter( alpha=0.8, size=2)+
theme_bw() +
scale_color_gradientn(colours = pal) +
scale_x_continuous(limits = c(0.6, 1.4)) +
scale_y_continuous(limits = c(0.6, 1.4)) +
geom_abline(slope = 1, col="gray50",linetype = 2)+
labs(x="RR for males", y="RR for females") +
theme(legend.position="none")
b<-COVID %>%
dplyr::select(SA3_CODE16,Population,RR,IRSAD) %>%
filter(Population=="Age 12-17" | Population=="Age 18-25") %>%
pivot_wider(names_from=Population,values_from=RR) %>%
mutate(compare=`Age 18-25` >`Age 12-17`) %>%
ggplot(aes(x=`Age 18-25`, y=`Age 12-17`, colour=IRSAD)) +
geom_jitter( alpha=0.8, size=2)+
theme_bw() +
scale_color_gradientn(colours = pal) +
scale_x_continuous(limits = c(0.6, 1.4)) +
scale_y_continuous(limits = c(0.6, 1.4)) +
geom_abline(slope = 1, col="gray50",linetype = 2)+
labs(x="RR for age 12-17", y="RR for age 18-25")+
theme(legend.position="none")
g_legend<-function(a.gplot){
tmp <- ggplot_gtable(ggplot_build(a.gplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
legend<-a +
theme(legend.position="bottom") +
theme(legend.key.height= unit(0.5, 'cm'),
legend.key.width= unit(1.2, 'cm'))
legend <- g_legend(legend)
grid.arrange(ggarrange(a,b,ncol=2) ,
legend, nrow=2,heights=c(10, 1))
Pooled effect
res <- results %>%
filter(Terms=="COVID") %>%
group_by(Population) %>%
group_modify(~ broom::tidy(extract_meta(data = .x))) %>%
pivot_wider(names_from=names,values_from=x) %>%
mutate('RR (95% CI)' = paste0(format(round(RR, digits=3), nsmall=3), " (" ,
format(round( CIL, digits = 3),nsmall=3)," - " ,
format(round(CIH, digits = 3),nsmall=3),")"))
res$y_loc=5:1
#main plot in the middle
plot1<- ggplot(res, aes(x = RR, y = y_loc, colour ) )+
geom_point() +
geom_errorbarh(aes(xmin = CIL, xmax = CIH), height = 0.3) +
ylab(NULL) +
scale_y_continuous(limits = c(0.8,5.5)) +
ggtitle(" ") +
xlab("Estimated RR (95% CI)") +
scale_colour_manual(values = colours ) +
theme(plot.title = element_text(hjust =1 ,size = 12,face = "bold"),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.border = element_blank(),
axis.line.x = element_line(colour = "black"),
panel.background = element_blank(),
panel.grid = element_blank(),
plot.margin = margin(0, 0, 0, 0, "cm"),
legend.position = "none") +
geom_vline(aes(xintercept = 1), linetype = "dashed", col="gray50")
#base plot for adding text
tab_base <- ggplot(res, aes(y = y_loc)) +
ylab(NULL) + xlab("") +
scale_y_continuous(limits = c(0.8,5.5)) +
theme(plot.title = element_text( size = 12), ## centering title on text
axis.text.x = element_text(color = "white"), ## need text to be printed so it stays aligned with figure but white so it's invisible
axis.line = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.y = element_blank(),
plot.margin = margin(0, 0, 0.3, 0, "cm"),
legend.position = "none",
panel.background = element_blank(),
panel.grid = element_blank()) +
scale_colour_manual(values = colours ) +
labs(x=" ")
# Predictors
tab1 <- tab_base +
geom_text(aes(x = 0, label = Population, y = y_loc), hjust = 0,nudge_x = 0.05) +
scale_x_continuous(limits = c(0,0.4)) +
ggtitle("Outcomes") +
theme(plot.title = element_text(hjust = 0.3))
# OR
OR1 <- tab_base +
geom_text(aes(x = 0, label = `RR (95% CI)`) , hjust = 0.5) +
ggtitle(" OR (95% CI)")+
theme(plot.title = element_text(hjust = 0.5))
gridExtra::grid.arrange(tab1, plot1, OR1,
layout_matrix = matrix(c(rep(1,2),
rep(2,6),
rep(3,3)), nrow = 1))
Meta regression
res <- results %>%
filter(Terms=="COVID") %>%
mutate(IRSAD= (IRSAD-995.5)/100) %>%
group_by(Population) %>%
group_modify(~reg_meta(data = .x)) %>%
mutate('RR ' = format(round(RR, digits=3), nsmall=3),
"95%CI" =paste0(format(round( CIL, digits = 3),nsmall=3)," - " ,
format(round(CIH, digits = 3),nsmall=3)),
'p-value'=ifelse(p<=0.001, "<0.001", format(round(p, digits=3), nsmall=3))) %>%
dplyr::select(-RR,-CIL,-CIH,-p) %>%
pivot_longer(names_to="Estimate", values_to="Value",-c("Population","Variable") ) %>%
pivot_wider(names_from = Population,values_from=Value)## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## NULL
## debug at <text>#7: print(summary(a))
##
## Mixed-Effects Model (k = 327; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 247.4643 -494.9285 -484.9285 -466.0403 -484.7393
##
## tau^2 (estimated amount of residual heterogeneity): 0.0074 (SE = 0.0009)
## tau (square root of estimated tau^2 value): 0.0861
## I^2 (residual heterogeneity / unaccounted variability): 68.96%
## H^2 (unaccounted variability / sampling variability): 3.22
## R^2 (amount of heterogeneity accounted for): 22.36%
##
## Test for Residual Heterogeneity:
## QE(df = 323) = 1020.5763, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 323) = 21.3740, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.0316 0.0092 3.4470 0.0006 0.0136 0.0497 ***
## IRSAD 0.0786 0.0106 7.4014 <.0001 0.0577 0.0995 ***
## StateVictoria -0.0254 0.0147 -1.7243 0.0856 -0.0543 0.0036 .
## RegionRegional 0.0220 0.0146 1.5132 0.1312 -0.0066 0.0507
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## NULL
## debug at <text>#7: print(summary(a))
##
## Mixed-Effects Model (k = 328; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 348.7922 -697.5843 -687.5843 -668.6806 -687.3957
##
## tau^2 (estimated amount of residual heterogeneity): 0.0032 (SE = 0.0004)
## tau (square root of estimated tau^2 value): 0.0569
## I^2 (residual heterogeneity / unaccounted variability): 61.52%
## H^2 (unaccounted variability / sampling variability): 2.60
## R^2 (amount of heterogeneity accounted for): 43.50%
##
## Test for Residual Heterogeneity:
## QE(df = 324) = 847.5023, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 324) = 40.4836, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.0378 0.0066 5.7557 <.0001 0.0249 0.0508 ***
## IRSAD 0.0503 0.0076 6.6052 <.0001 0.0353 0.0653 ***
## StateVictoria 0.0966 0.0112 8.6598 <.0001 0.0747 0.1185 ***
## RegionRegional 0.0315 0.0105 3.0002 0.0029 0.0109 0.0522 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## NULL
## debug at <text>#7: print(summary(a))
##
## Mixed-Effects Model (k = 327; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 349.4359 -698.8718 -688.8718 -669.9836 -688.6826
##
## tau^2 (estimated amount of residual heterogeneity): 0.0030 (SE = 0.0004)
## tau (square root of estimated tau^2 value): 0.0545
## I^2 (residual heterogeneity / unaccounted variability): 55.25%
## H^2 (unaccounted variability / sampling variability): 2.23
## R^2 (amount of heterogeneity accounted for): 42.91%
##
## Test for Residual Heterogeneity:
## QE(df = 323) = 732.5255, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 323) = 38.8790, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.0838 0.0065 12.8567 <.0001 0.0710 0.0966 ***
## IRSAD 0.0665 0.0076 8.6920 <.0001 0.0514 0.0815 ***
## StateVictoria 0.0655 0.0111 5.8995 <.0001 0.0437 0.0874 ***
## RegionRegional 0.0323 0.0105 3.0947 0.0021 0.0118 0.0529 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## NULL
## debug at <text>#7: print(summary(a))
##
## Mixed-Effects Model (k = 328; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 294.2438 -588.4875 -578.4875 -559.5838 -578.2988
##
## tau^2 (estimated amount of residual heterogeneity): 0.0046 (SE = 0.0006)
## tau (square root of estimated tau^2 value): 0.0682
## I^2 (residual heterogeneity / unaccounted variability): 65.46%
## H^2 (unaccounted variability / sampling variability): 2.89
## R^2 (amount of heterogeneity accounted for): 12.17%
##
## Test for Residual Heterogeneity:
## QE(df = 324) = 923.1324, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 324) = 10.1745, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt -0.0497 0.0079 -6.3202 <.0001 -0.0652 -0.0343 ***
## IRSAD 0.0457 0.0090 5.0922 <.0001 0.0280 0.0633 ***
## StateVictoria 0.0148 0.0126 1.1763 0.2403 -0.0099 0.0395
## RegionRegional 0.0155 0.0125 1.2427 0.2149 -0.0090 0.0401
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)
## Called from: reg_meta(data = .x)
## debug at <text>#6: print(unique(data$population))
## NULL
## debug at <text>#7: print(summary(a))
##
## Mixed-Effects Model (k = 331; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 380.8468 -761.6936 -751.6936 -732.7438 -751.5067
##
## tau^2 (estimated amount of residual heterogeneity): 0.0028 (SE = 0.0004)
## tau (square root of estimated tau^2 value): 0.0533
## I^2 (residual heterogeneity / unaccounted variability): 63.02%
## H^2 (unaccounted variability / sampling variability): 2.70
## R^2 (amount of heterogeneity accounted for): 37.81%
##
## Test for Residual Heterogeneity:
## QE(df = 327) = 876.0617, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 327) = 34.7422, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.0372 0.0060 6.1879 <.0001 0.0254 0.0490 ***
## IRSAD 0.0592 0.0070 8.5049 <.0001 0.0455 0.0729 ***
## StateVictoria 0.0505 0.0101 5.0069 <.0001 0.0306 0.0703 ***
## RegionRegional 0.0256 0.0096 2.6728 0.0079 0.0068 0.0445 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## debug at <text>#8: ret <- data.frame(cbind(exp(cbind(a$b, a$ci.lb, a$ci.ub)), a$pval))
## debug at <text>#9: names(ret) <- c("RR", "CIL", "CIH", "p")
## debug at <text>#10: ret$Variable <- c("Intercept", "IRSAD", "Victoria", "Regional")
## debug at <text>#11: return(ret)res$Variable[c(2:3,5:6,8:9,11:12)]<-""
print_table(res, align = c("l", "l",rep("c",5)), colwidth = c(2,3,rep(4.4,5)))## <caption> (\#tab:)</caption>Variable | Estimate | Age 12-17 | Age 18-25 | Females | Males | All young people |
Intercept | RR | 1.032 | 1.039 | 1.087 | 0.951 | 1.038 |
95%CI | 1.014 - 1.051 | 1.025 - 1.052 | 1.074 - 1.101 | 0.937 - 0.966 | 1.026 - 1.050 | |
p-value | <0.001 | <0.001 | <0.001 | <0.001 | <0.001 | |
IRSAD | RR | 1.082 | 1.052 | 1.069 | 1.047 | 1.061 |
95%CI | 1.059 - 1.105 | 1.036 - 1.067 | 1.053 - 1.085 | 1.028 - 1.065 | 1.047 - 1.076 | |
p-value | <0.001 | <0.001 | <0.001 | <0.001 | <0.001 | |
Victoria | RR | 0.975 | 1.101 | 1.068 | 1.015 | 1.052 |
95%CI | 0.947 - 1.004 | 1.078 - 1.126 | 1.045 - 1.091 | 0.990 - 1.040 | 1.031 - 1.073 | |
p-value | 0.086 | <0.001 | <0.001 | 0.240 | <0.001 | |
Regional | RR | 1.022 | 1.032 | 1.033 | 1.016 | 1.026 |
95%CI | 0.993 - 1.052 | 1.011 - 1.054 | 1.012 - 1.054 | 0.991 - 1.041 | 1.007 - 1.046 | |
p-value | 0.131 | 0.003 | 0.002 | 0.215 | 0.008 |